Finite-size effects in the interfacial stiffness of rough elastic contacts 
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The stiffness K of the interface between two rough but nominally flat elastic solids is studied 
analytically and numerically. Over a wide range of applied squeezing pressure p, K rises linearly 
with p. Sublinear power law scaling is observed at small p, but computer simulations are used to 
demonstrate that this is a finite-size effect. Moreover, we derive accurate, analytical expressions for 
the exponents and prefactors of the low-pressure scaling of K by extending the contact mechanics 
theory of Persson to systems of finite size. 



Two solids in mechanical contact tend to touch at only 
a miniscule fraction of their apparent contact area Aq, be- 
cause their surfaces are microscopically rough [2-IH . This 
imperfect contact has profound implications for transmis- 
sion of charge, heat and forces through the interface. The 
effect of the interface can be expressed in terms of an in- 
terfacial stiffness or conductance that adds in series with 
the bulk response of two solids with ideal fiat surfaces 
[5], IH . Improved theories of these interfacial contributions 
are important because they frequently dominate the to- 
tal response of the system and are a strong function of 
the normal force F (or load) pushing the solids together. 
In this paper we consider the scaling of stiffness with F 
for nonadhesive self-affine rough surfaces. The results 
are more generally applicable since the shear and normal 
stiffness and electrical and heat conductance are all pro- 
portional to each other within linear-response continuum 
mechanics [6j • 

In a pioneering experimental work, Berthoud and 
Baumberger found that the interfacial stiffness was pro- 
portional to F for nonadhesive solids with very differ- 
ent elastic properties The proportionality can be ex- 
pressed as 



K = p/u , 



(1) 



where K is the interfacial stiffness normalized by Aq , and 
p = F/Aq. The characteristic length uq was found to be 
of order the combined root mean squared (rms) rough- 
ness h rms (~ 1/im) of the surfaces. The surfaces had 
sclf-affinc fractal roughness that is common in experi- 
ments. Berthoud and Baumberger rationalized their ob- 
servations within the contact mechanics theory of Green- 
wood and Williamson [7], which, however, is based on 
hypotheses that later turned out to be unjustified 
Nonetheless, the results of additional experiments 111 , 
and computer simulations of elastic contacts (TlL Tl3- 
with self-affine, fractal roughness, are consistent with 
Eq. ((T|). Moreover, the proportionality coefficient uq 
agrees, to within 0(10%), with uq « 0.4/i rms , derived 
from the parameter-free contact mechanics theory of 
Persson fia[l7|. 



The interfacial stiffness can be determined from the to- 
tal stiffness K to t and the stiffnesses of ideal fiat bounding 
solids K\ and K% using the rule for springs in series 



K- l =K^-K?-K? 



(2) 



An alternative approach is to measure — or to compute 
— the mean interfacial separation u. Changes in u are 
a direct measure of the deformation attributable to the 
interface and K = -dp/du where the sign reflects the fact 
that u decreases with increasing confining force. In the 
range of validity of ([1]), this differential relation can be 
solved to yield another testable prediction 



p = p exp(-u/u ), 



(3) 



where po is an integration constant. Persson theory finds 
that po = (3E* , where E* is the effective clastic modulus 
and (3 is dimensionless. Like uq, ft only depends on the 
spectral properties of the surface [8]. Analytical expres- 
sions for Mo and /3 and computer simulations agree again 
to within 0(10%) 

In a recent Letter [20j|7pohrt and Popov challenged the 
established results on interfacial stiffness by proposing a 
sublinear K oc p a power law deduced from numerical sim- 
ulations. Specifically, they reported a = 0.2567 x (3 -if), 
where H is the Hurst roughness exponent. This estimate 
was later corrected to a = 0.266 x (3 - H) and scaling 
arguments were presented for a = 1/(1 + H) [2l|. Pohrt 
et al. argued that their results differed from previous 
ones because their surfaces were "truly fractal' [2lJ, i.e., 
roughness lived on wavelengths all the way to the linear 
system size L. In particular they state: "Whenever the 
surfaces are truly fractal with no cut-off wavelength, a 
power law applies" [2lj |. 

In this Letter, we unravel the origin of the discrepancy 
between the established results and the new findings. To 
do so, we analyze finite-size effects in numerical simula- 
tions. In addition, we derive analytical expressions, free 
of adjustable parameters, to estimate finite-size effects. 
While the essence of the calculations is presented in the 
main part of this work, details, such as the derivation 



of prefactors, can be found in the supplementary ma- 
terial [22J. It also contains unpublished experiments in 
support of Persson's contact mechanics theory. 

We first summarize the arguments for how Eq. (fTJ) 
arises from the self-affinity of interfaces Q. The key 
idea is that when there are a large number of separated 
contacting patches, the distribution of contacts is self- 
similar. As the load increases, existing contact patches 
grow and new, small contacts are formed. This hap- 
pens in such a way that the distributions of contact sizes 
and local pressures remain approximately constant over 
a wide range of loads 0, [23j . An immediate consequence 
is a linear relation between real contact area A and p, 
which has been confirmed in many simulations, includ- 
ing all numerical studies cited here. The spatial correla- 
tions between contacting areas and local stresses are also 
the same up to a prefactor that grows linearly with load 
because of a sum rule 2J]. Since the system responds 
linearly, the elastic energy U e \ is given by an integral of 
an elastic Greens function times the Fourier transform 
of the stress-stress correlation function and must thus be 
proportional to load: 



U c \ = u A p. 



(4) 



Since the elastic energy is equal to the work done by the 
external load (assuming hard-wall interactions and no 
adhesion), it follows that dU c \ = u^A^dp = -Aop(u)du. 
This last relation is identical to (J3|) and thus also to ((T|). 

When p is so small that two finite surfaces start touch- 
ing, the interface cannot yet behave in a self-similar fash- 
ion. The reason is that contact occurs only near the 
highest asperity whose height determines the separation 
at first contact u c . As a consequence, the validity of the 
arguments leading to @ and thus to ([1} — or any the- 
ory valid in the thermodynamic limit — cannot hold at 
small p. As already pointed out earlier, finite-size effects 
then become important [25| . Specifically, for a finite sys- 
tem p vanishes for (finite) u > u c , while for an infinite 
system p is always non- vanishing. Thus, p must initially 
decay faster with increasing u in a finite system than in 
an infinite system where Eq. (TTJ) holds. In the opposite 
case of large p, a finite system approaches complete con- 
tact, u = 0, at finite pressure but infinite systems do not 
because they have infinitely deep valleys. Thus, the in- 
terfacial stiffness of small systems at large p must diverge 
more quickly than that of large ones. One may conclude 
that contact formation of the highest peak and the low- 
est valley depend on the specific realization of a surface. 
However, for intermediate pressures, universal behavior 
may be found as long as the roughness has well-defined 
statistical properties. 

To study finite-size effects, we performed large-scale 
numerical simulations of nonadhesive contact between 
a rigid self-affine surface and an isotropic elastic sub- 
strate with effective modulus E* and Poisson number 
0.5 using well-established methods [2(| 27 j. Surfaces 




FIG. 1: (Color online) Log-log plot of the nondimensional 
contact stiffness Kh Ilns /E* vs. nondimensional pressure p/E* 
for self-affine fractal surfaces with H = 0.7 and rms slope 
0.1. In all cases the surface is resolved with 8192 points in 
each direction, L/X\ = 4096 and the ratio of system size to 
roll-off wavelength, L/X r is indicated. The (red) solid line 
is the prediction of Persson's theory while the dashed (red) 
line is the linear regime. Inset: Nondimensional pressure vs 
separation for the same surfaces. 



were self-affine with Hurst exponent H between a short 
wavelength cut-off Ai and long-wavelength roll-off A r (see 
Fig. 1 in [22J). The amplitudes of the Fourier transforms 
for the height /i(q) were drawn from a Gaussian distri- 
bution. Their variances reflect the roughness spectrum 
C(q) for each reciprocal space vector q. For wavenum- 
bers qo < q < q r (where go = w/L and q r = 7r/A r ), we 
used C(q) = C(q r ), while the spectrum reflected the de- 



-2(1+H) 



the 



sired self-affine scaling C(q) = C(q r )(q/q r )~ 
range q r < q < q\ (where qi = n/Xi). Periodic boundary 
conditions with period L were applied in the plane of the 
substrate. 

Figure [T] shows typical results for the contact stiffness 
versus pressure. There is a linear relation at intermediate 
loads in all cases and a more rapid rise of K with p as full 
contact is approached. Both regimes are well-described 
by Persson's contact mechanics theory (red line), which 
requires only the surface roughness power spectrum and 
the effective modulus as input. We also find a transition 
to power law scaling at low loads. The prefactor to the 
power law varies by as much as an order of magnitude 
from one specific realization to the next (not shown ex- 
plicitly). The transition from power law to a linear K(p) 
behavior is particularly sensitive to the magnitude of a 
few random Fourier components at the smallest wavevec- 
tors as well as to their relative phases. The value of u c 
is also very sensitive to these Fourier components and 
decreases with L/X r . It cuts off the exponential relation 
between p and K shown in the inset. 

Even for the case where \ r /L = 1, the results in Fig- 
ure [T] follow linear scaling (Eq. ([T])) for more than one 
decade. The range of validity rapidly extends to lower p 
as L/X r increases. Thus, the more closely the thermody- 



namic limit is approached — or the more significant the 
statistical distribution of contacting peaks — the more 
accurate is Eq. (J]). Given typical A r , e.g., O(10^m) 
for polished steel and 0(1 cm) for asphalt, one can see 
that power law scaling matters only for either very small 
systems or very small loads. 

In this low load regime, Fig. Q] indicates K oc p a be- 
havior with a « 0.6 for H = 0.7. Thus, our sm all-p ressure 
results for L/\ r = 1 are consistent with Refs. 20j,|21|. In 
the following, we propose an explanation for this power 
law by incorporating the estimation of finite-size effects 
into Persson's contact mechanics theory. The goal is 
to find an expression for the elastic energy because it 
allows us to calculate the contact stiffness. We reex- 
press a small change of the elastic energy dU c \ = -pA du 
as dU e i = -pAodp(du/dp). Inserting K = -dp/du and 
F = pA yields 



p = K- 



-dUei 
dF 



(5) 



Our approach is motivated by the fact that the elas- 
tic energy is dominated by the longest wavelength modes 
For a single contacting region around the highest 
peak, the longest wavelength will scale with the radius ro 
of the smallest circle that encloses the contacts. We will 
first calculate the elastic energy for a single Hertzian-like 
mesoscale asperity with radius of curvature R and con- 
tact radius ro. Then we show that including roughness 
at smaller wavelengths gives the same power law scaling. 

An effective asperity radius is calculated from the 
roughness at scales larger than ro- The local curvature 
V 2 h corresponds to q 2 h(q) in Fourier space. Thus R can 
be estimated as: 

— oc / d 2 q\q 2 h(q)\ 2 oc / dq q 5 C(q). (6) 

For self-affine fractal roughness the surface roughness 
power spectrum C(q) oc q- 2 ( 1+H ) . This gives R oc r]f H , 
where we have assumed that the lower integration bound 
to the last integral must be negligible at a small load. 
This condition is fulfilled as long as ro « A r . 

According to Hertzian contact mechanics, ro oc 
(RF) 1 / 3 . Inserting R oc r^ H and solving for ro, we ob- 
tain 



r o « 



The elastic energy stored within a Hertzian contact is 
oc F8 where the penetration depth i5 oc r 2 /R oc Tq. 
We obtain 



[7 (0) oc ^(i+2H)/(i+H) 



and from Eq. ([5]) 



(7) 



(8) 



We now show that the elastic energy U^p due to mi- 
croscale roughness within the mesoscale asperity also 
scales with F^ 1+2H ^^ 1+H \ The main assumption now is 
that the contact pressure within the mesoscale asperity 
contact region is high enough that the contact mechan- 
ics theory by Persson can be applied. Then from Q, 
= uiAipi, where A\ = irr 2 . is the (nominal) contact 
area at the mesoscale and p\ = F\A\. The term u\ is of 
order the rms roughness including only roughness com- 
ponents with wavelength A < ro. This can be written 
as 



h 2 . ms = 2n dq qC(q) oc (*) ™ - (f ) 

Since Ai « ro (unless H is close to 0) one obtains h rms oc 
and ui oc Inserting r oc F 1 ^ 1+H \ we get U j oc 

F a + 2H)/(i+H) as in Eq ^ Q 

From the above treatment we predict that the stiffness 
K scales as p a with a = 1/(1 + H). Fig. [5] shows K(p) 
relations obtained numerically in the finite-size regime 
for different values of H. Rough estimates for a were 
obtained by fitting to the lowest four data points. The 
results from the simulations are: a(H = 0.3) = 0.72 (see 
also below), a(H = 0.5) = 0.66, and a(H = 0.7) = 0.59. 
These values compare well to the theoretical predictions, 
a(H = 0.3) * 0.769, a(H = 0.5) « 0.667, and a(H = 0.7) « 
0.588, particularly if one keeps in mind that calculations 
become more approximate as H approaches zero. 
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FIG. 2: (Color online) Dimensionless interfacial stiffness 
Kh Ilna /E* as a function of pressure p/E* in the finite-size 
region for L/Ai = 4096 and different Hurst exponents H. All 
calculations are for L/\ r = 1. Solid lines show a fit to the first 
four data points. 

We note here that the expression for the microasper- 
ity contribution to the total elastic energy depends on 
the elastic coupling between the asperities. Any deriva- 
tion neglecting this coupling [2(| cannot describe the cor- 
rect physics, even if the resulting scaling is similar to 
K oc p l li 1+H ) . Moreover, probing the constitutive rela- 
tion between pressure and stiffness at a mesoscale will 
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FIG. 3: (Color online) Exponent a as a function of 



(Ai/A r ) u 



Theoretically predicted value is denoted by 



The inset shows selected numerical data for K(p) from which 
a was deduced. Data for the following values of A,-/Ai is pre- 
sented: 4096 (•), 2048 (♦), 768 (*), 512 (A). 



entail much larger fluctuations than in a multi-asperity 
contact at the same pressure but larger value of L/X r . 
The arguments that lead to K oc p l IK lJtH ) hold when 



Ai « r « A. r . Since R oc 7q , the radius of the 
mesoasperity diverges as the contact area grows and 
r o <V- br this limit, the mesoasperities are flat, both 
Eqns. ([U and (O hold, and we rediscover the thermo- 
dynamic limit K oc p 13- 19j. On the other hand, if 
ro < Ai the surface of the mesoasperity is smooth. The 
upper integration bound in (|6|) is then given by the short 
wavelength cutoff q\ = n/X\ and R is constant. This ul- 
timately must lead to traditional Hertz behavior where 
K oc p 1 / 3 . Indeed, we find that reducing the ratio of 
A r /Ai gives an exponent that approaches the one ex- 
pected from Hertz contact mechanics. This additional 
scaling regime limits the range where a should be ob- 
served and complicates its measurement in simulations. 
Fig. [3] shows the results of an attempted extrapolation 
to the "fractal limit" Ai/A r -> for the value of H = 0.3, 
which had the largest discrepancy between theory and 
simulation. Despite quite large stochastic scatter, we 
conclude that the value of a = 1/(1 + H) is consistent 
with the simulations. 

Finally, we address how the finite size power-law region 
depends on linear system size L and roll-off length A r . 
Following along the lines of the above derivation, it is 
straightforward to compute the full expression for the 
interfacial stiffness: 



K(p) = 9 



E X r 



pL 2 q r \ 
E*h^ s s l l 2 } 



1/(1+H) 



(9) 



where 1/s = 1 + H[l - (qo/q r ) 2 ]- The prefactor 8 depends 
only on the Hurst exponent H , but for H > 0.3 variation 
is restricted to 1.0 < 9 < 1.3. The full details of this calcu- 



lation can be found in the supplementary material [22j . 
By equating ([9]) with K(p) = p/-fh lms , where 7 « 0.4, 
we obtain an estimate for the pressure p c at which the 
stiffness crosses over from power-law to linear behavior: 



E 



. q r h TI 



-1/2H 



(07) 



(1+H)/H 



(10) 



For different realizations of the surface we expect large 
fluctuations in the prefactor of the power-law relation ((9]) , 
and hence also in the cross-over pressure p c . Neverthe- 
less, for the data shown in Fig. [T] we find p c /E* « 8x 10~ 5 
for q r /qo = 1 and p c /E* « 4x 10~ 6 for q r /qo = 8, in ex- 
cellent agreement with the numerical data. Generally, 
the cross-over pressure p c decreases with increasing linear 
system size L. Equation (|10p also reveals the importance 
of separation between L and the roll-off length A,- . Scale 
separation pushes the crossover to lower pressure even 
more rapidly since the ratio L/X r enters quadratically. 
In the thermodynamic limit L/X r -*■ 00, the power-law 
region vanishes alltogether. 

We conclude that the previously reported K oc p and 
p oc exp(-u/uo) laws are satisfied when there is a statis- 
tical ensemble of large peaks in contact. While there is a 
region where K scales sublinearly with p, prefactor and 
onset, u c , and even the slope parameter a, are not uni- 
versal and have no thermodynamic limit. Instead, they 
have large fluctuations from one sample to the next, even 
though surface topographies are constructed using iden- 
tical procedures. 
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